%% Spectrum Math

close all
clear all
clc

a = csvread('illum_a_spec.csv');
d50 = csvread('illum_d50_spec.csv');
d65 = csvread('illum_d65_spec.csv');
tl84 = csvread('illum_tl84_spec.csv');

ave_a = sum(a(:,2:end),2)/5;
ave_d50 = sum(d50(:,2:end),2)/5;
ave_d65 = sum(d65(:,2:end),2)/5;
ave_tl84 = sum(tl84(:,2:end),2)/5;

norm_a = ave_a/sum(ave_a);
norm_d50 = ave_d50/sum(ave_d50);
norm_d65 = ave_d65/sum(ave_d65);
norm_tl84 = ave_tl84/sum(ave_tl84);

ymax = max([max(norm_a);max(norm_d50); max(norm_d65);max(norm_tl84)]);

wavelength_a = a(:,1);
wavelength_d50 = d50(:,1);
wavelength_d65 = d65(:,1);
wavelength_tl84 = tl84(:,1);

%% Plot Spectral Data
close all

figure
hold on
plot(wavelength_a, norm_a, 'r', 'LineWidth', 2);
plot(wavelength_d50, norm_d50, '-b', 'LineWidth', 2);
plot(wavelength_d65, norm_d65, '--b', 'LineWidth', 2);
plot(wavelength_tl84, norm_tl84, 'k', 'LineWidth', 2);
title('Spectral Power Distribution of Sources');
xlabel('Wavelength (nm)');
ylabel('Area Normalized Radiance');
axis([350 800 0 ymax]);
legend('A', 'D50', 'D65', 'TL84', 'Location', 'NorthEast');
hold off

print( '-depsc','spec_all_source.eps');

%% Plot Comparison (Fluorescent)

figure

hold on
plot(wavelength_d50, norm_d50, '-b', 'LineWidth', 2);
plot(wavelength_d65, norm_d65, '--b', 'LineWidth', 2);
title('Spectral Power Distribution of Two Fluorescent Sources');
xlabel('Wavelength (nm)');
ylabel('Area Normalized Radiance');
axis([350 800 0 ymax]);
legend('D50', 'D65', 'Location', 'NorthEast');
hold off

print( '-depsc','spec_d50_d65.eps');

%% Plot Comparison (A and Fluorescent)

figure

hold on
plot(wavelength_a, norm_a, 'r', 'LineWidth', 2);
plot(wavelength_d65, norm_tl84, 'k', 'LineWidth', 2);
title('Spectral Power Distribution of a Blackbody Radiator and a Fluorescent Ligh Source');
xlabel('Wavelength (nm)');
ylabel('Area Normalized Radiance');
axis([350 800 0 ymax]);
legend('A', 'D65', 'Location', 'NorthEast');
hold off

print( '-depsc','spec_a_tl84.eps');

%% Color Plotting Math

load('color_matching_func.mat');

cmf = color_matching_func(21:5:421,:);

xbar = cmf(:,2);
ybar = cmf(:,3);
zbar = cmf(:,4);

% Illuminant A
X_a = sum(xbar.*norm_a);
Y_a = sum(ybar.*norm_a);
Z_a = sum(zbar.*norm_a);

x_a = X_a/(X_a + Y_a + Z_a);
y_a = Y_a/(X_a + Y_a + Z_a);

% Illuminant D50
X_d50 = sum(xbar.*norm_d50);
Y_d50 = sum(ybar.*norm_d50);
Z_d50 = sum(zbar.*norm_d50);

x_d50 = X_d50/(X_d50 + Y_d50 + Z_d50);
y_d50 = Y_d50/(X_d50 + Y_d50 + Z_d50);

% Illuminant D65
X_d65 = sum(xbar.*norm_d65);
Y_d65 = sum(ybar.*norm_d65);
Z_d65 = sum(zbar.*norm_d65);

x_d65 = X_d65/(X_d65 + Y_d65 + Z_d65);
y_d65 = Y_d65/(X_d65 + Y_d65 + Z_d65);

% Illuminant TL84
X_tl84 = sum(xbar.*norm_tl84);
Y_tl84 = sum(ybar.*norm_tl84);
Z_tl84 = sum(zbar.*norm_tl84);

x_tl84 = X_tl84/(X_tl84 + Y_tl84 + Z_tl84);
y_tl84 = Y_tl84/(X_tl84 + Y_tl84 + Z_tl84);

% Spectrum Locus

xlocus = xbar./(xbar + ybar + zbar);
ylocus = ybar./(xbar + ybar + zbar);

% Black Body Curve

bbfirst = [2000:7000]';
[rowfirst,~] = size(bbfirst);
xDfirst = -4.6070*(repmat(10^9,rowfirst,1)./(bbfirst.^3)) + 2.9678*(repmat(10^6,rowfirst,1)./(bbfirst.^2)) + 0.09911*(repmat(10^3,rowfirst,1)./bbfirst) + 0.244063;
bbsecond = [7000:25000];
[rowsecond,~] = size(bbsecond);
xDsecond = -2.0064*(repmat(10^9,rowsecond,1)./(bbsecond.^3)) + 1.9018*(repmat(10^6,rowsecond,1)./(bbsecond.^2)) + 0.24748*(repmat(10^3,rowsecond,1)./bbsecond) + 0.237040;

yDfirst = -3*xDfirst.^2 + 2.87*xDfirst - 0.275;
yDsecond = -3*xDsecond.^2 + 2.87*xDsecond - 0.275;

%% Plot Spectral Locus and xy coordinates

figure

hold on
plot(x_a, y_a, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
plot(x_d50, y_d50, 'o', 'MarkerFaceColor', 'g', 'MarkerSize', 8);
plot(x_d65, y_d65, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 8);
plot(x_tl84, y_tl84, 'o', 'MarkerFaceColor', 'k', 'MarkerSize', 8);
plot(xlocus, ylocus, 'k', 'LineWidth', 2);
plot([xlocus(1),xlocus(end)], [ylocus(1), ylocus(end)], 'k', 'LineWidth', 2);
plot(xDfirst, yDfirst,'-k');
plot(xDsecond, yDsecond, '-k');
title('Chromaticity Plot of Sources');
xlabel('x');
ylabel('y');
axis([0 0.9 0 0.9]);
set(gca,'YTick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
set(gca,'XTick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
box on
axis square
legend('A', 'D50', 'D65', 'TL84', 'Location', 'NorthEast');
hold off

print( '-depsc','xy_chromaticity.eps');

%% u'v' Calculation

uprime_a = (4*X_a)/(X_a + 15*Y_a + 3*Z_a);
vprime_a = (9*Y_a)/(X_a + 15*Y_a + 3*Z_a);

% Illuminant D50

uprime_d50 = (4*X_d50)/(X_d50 + 15*Y_d50 + 3*Z_d50);
vprime_d50 = (9*Y_d50)/(X_d50 + 15*Y_d50 + 3*Z_d50);

% Illuminant D65

uprime_d65 = (4*X_d65)/(X_d65 + 15*Y_d65 + 3*Z_d65);
vprime_d65 = (9*Y_d65)/(X_d65 + 15*Y_d65 + 3*Z_d65);

% Illuminant TL84

uprime_tl84 = (4*X_tl84)/(X_tl84 + 15*Y_tl84 + 3*Z_tl84);
vprime_tl84 = (9*Y_tl84)/(X_tl84 + 15*Y_tl84 + 3*Z_tl84);

% Spectrum Locus

uprimelocus = (4*xbar)./(xbar + 15*ybar + 3*zbar);
vprimelocus = (9*ybar)./(xbar + 15*ybar + 3*zbar);

% Blackbody Curve Calculation HERE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot u'v' chromaticity plot

figure

hold on
plot(uprime_a, vprime_a, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
plot(uprime_d50, vprime_d50, 'o', 'MarkerFaceColor', 'g', 'MarkerSize', 8);
plot(uprime_d65, vprime_d65, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 8);
plot(uprime_tl84, vprime_tl84, 'o', 'MarkerFaceColor', 'k', 'MarkerSize', 8);
plot(uprimelocus, vprimelocus, 'k', 'LineWidth', 2);
plot([uprimelocus(1),uprimelocus(end)], [vprimelocus(1), vprimelocus(end)], 'k', 'LineWidth', 2);
% plot(xDfirst, yDfirst,'-k');
% plot(xDsecond, yDsecond, '-k');
title('Chromaticity Plot of Sources');
xlabel('u''');
ylabel('v''');
axis([0 0.7 0 0.7]);
set(gca,'YTick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7])
set(gca,'XTick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7])
box on
axis square
legend('A', 'D50', 'D65', 'TL84', 'Location', 'SouthEast');
hold off

print( '-depsc','uprime_vprime_chromaticity.eps');